On the numerical solution of diffusion systems with localized, gradient-driven, moving sources

formulation Let us, before proceeding, for convenience first define the index sets of the dynamic fields, Jd = {1, . . . ,Md}, static (dynamic) fields, Js = {Md + 1, . . . ,M}, dynamic states, Id = {1, . . . , Nd}, and static (dynamic) states, Is = {Nd + 1, . . . , N}. The equations (6.2)–(6.6) together constitute the dynamical system behind the model that can be written in the form ż = f(t, z,y) 0 = g(t, z,y) with z(0) = z0 y(0) = y0. (6.8) Here, the z is composed of the dynamic fields and states, i.e., the fields ρj for j ∈ Jd and the states ui for i ∈ Id, and therefore consists of a selection of components of the total model state x in (6.1). Likewise, the vector y is composed of the static fields and states, i.e., the fields ρj for j ∈ Js and the states ui for i ∈ Is, forming the remaining part of the model state x. The initial condition consisting of the vectors z0 and y0 has to be chosen in such a way that it obeys g(0, z0,y0) = 0. The dynamical system (6.8) can be turned into a ‘lower-dimensional’ system for z only, if we assume that for given values of z and t we can solve y uniquely from g(t, z,y) = 0. This yields y as a function of t and z, y = h(t, z), resulting in the system ż = f(t, z, h(t, z)) with z(0) = z0. (6.9) 6.2. SIMULATION FRAMEWORK 95 We will now examine equation g(t, z,y) = 0 in more detail. Here, the vector z is composed of the dynamic fields and states, i.e., ρj with j ∈ Jd and ui with i ∈ Id. Assuming that these and the time t are given, solving y from g(t, z,y) = 0 comes down to solving the static fields and states from 0 = Ljρj + ∑ i∈I=Id∪Is σji(si)TriS, j ∈ Js, 0 = Gi ( t,ui,ρ(ri), ∂xρ(ri), ∂yρ(ri) ) , i ∈ Is. (6.10) In the case where we use the general functions Gi we have to solve this full system, which is nonlinear and infinite-dimensional. Using Newton iteration is a possibility, but one that requires solving elliptic equations every iteration step. An appealing alternative exists if we use the extra assumption that the functions Gi are of the special form (6.7). We can then turn the system into a finite-dimensional system for which we do not have to solve elliptic equations for every iteration step, but we have to do it only once. To this end we first eliminate the elliptic equations from the system (6.10). Using these, the fields ρj can be expressed in the si and ri by writing ρj = − ∑ k∈Id∪Is σjk(sk)L−1 j TrkS, j ∈ Js =⇒ ρ(ri) = −diag ([σ(s̄)][S(r̄, ri)]) . (6.11) Here, the operators L−1 j , that commute with the scalars σjk, denote the inverse operators of Lj with respect to the boundary conditions (6.2). The s̄ and r̄ denote the vector (s1 , . . . , s T N ) T and (r1 , . . . , r T N ) T , respectively, while the matrices σ and S are defined by [σ(s̄)]jk = σjk(sk) and [S(r̄, ri)]kj = (L−1 j TrkS)(ri), respectively. The function diag(·) is defined to return the diagonal vector of its argument. Defining the matrix Sx by [Sx(r̄, ri)]kj = ∂x(L−1 j TrkS)(ri) yields a similar expression for ∂xρ(ri) with S replaced with Sx, while a similar definition of Sy results in a similar expression for ∂yρ(ri). The second equation in system (6.10) can now be written as 0 = Gi ( t, (ri, si),−diag ( [σ(s̄)][S(r̄, ri)] ) ,−diag ([σ(s̄)][Sx(r̄, ri)]), − diag ([σ(s̄)][Sy(r̄, ri)])), i ∈ Is, (6.12) where we wrote the state ui as (ri, si). In doing this we have replaced the infinitedimensional system (6.10) with the finite-dimensional system (6.12), where the s̄ and the r̄ are composed of all the si and ri, respectively (i ∈ Is ∪ Id), but only the static si, ri are the unknowns. Although the resulting system (6.12) is essentially finite-dimensional, applying Newton iteration requires solving elliptic equations each iteration step, needed for evaluation of the matrices S, Sx and Sy. However, if the Gi are of the form (6.7), the system decouples. It is then possible to solve the static ri 96 CHAPTER 6. NUMERICAL SOLUTION OF FRAMEWORK SYSTEMS first using the functions Gi . Afterwards, all ri are known and elliptic equations have to be solved to find the S, Sx and Sy. Finally, Newton iteration is used to solve 0 = Gi ( t, si,−diag ( [σ(s̄)][S] ) ,−diag ([σ(s̄)][Sx]),−diag ([σ(s̄)][Sy])), i ∈ Is, for the static si, where we left out the arguments of the matrices S, Sx and Sy. 6.3 Numerical methods in AGTools In this section we will go into the numerical machinery implemented in AGTools for approximating solutions of the systems (6.2)–(6.6). We will start with making some general remarks on the adopted approach for discretizing the model equations. In general there are two approaches for writing down full discretizations of time dependent PDE systems. The most used one is called the Method of Lines (MOL) approach and starts with a spatial discretization of the dependent fields and their differetial equations, turning the system in a large, but finite-dimensional ODE-system, called the semi-discrete system. Then a suitable time integrator is selected for the temporal discretization to yield the fully discrete solution. An advantage of the MOL approach is that one can choose a suitable method from a large collection of time integration methods for ODEs that are available. The downside is that the semi-discrete systems might become very complicated due to the presence of certain discretizationor interpolation operators. Direct application of, for example Rosenbrock methods, which involve the evaluation of Jacobians, can become cumbersome or even impossible. The second approach is the so-called Rothe approach [39]. Instead of first choosing a spatial discretization it starts by selecting a time integration method. This will result in a sequence of PDEs in time containing only spatial derivatives (boundary value problems). In this approach the PDEs are often stated as an abstract ODE in a certain Banach space making that the analysis of the used time integration methods moves to the realm of functional analysis and therefore becomes much more difficult The harder analysis however is accompanied by a number of advantages. First, the approach seems to have a cleaner appearance, not having to deal with difficult ODEs with discontinuities that are the result of spatial discretizations, but instead with elliptic equations coupled to algebraic equations, where everything is still smooth from spatial perspective. Second, it allows for nice error estimators, as is clearly described and illustrated by Lang [31]. We will adopt here the Rothe approach and not consider any functional analytic aspects, but take the practical approach in which we assume that our time integration methods work well for our cases, i.e., do not display any instability behavior. Time integration With respect to the time integration it is of importance that system (6.2)–(6.6) consists of stiff and non-stiff parts. The diffusion in the field equations gives rise to 6.3. NUMERICAL METHODS IN AGTOOLS 97 stiffness, while the nonlinear state equations are not necessarily stiff. For the time integration of stiff equations we would like to use an implicit method, but using such a method becomes very complicated for the system at hand and is not really suitable for the non-stiff part as well. IMEX scheme A class of schemes that seems to be appropriate here, is the class of IMEX (IMplicit/EXplicit) schemes. Such schemes can be applied to dynamical systems of the form ż = F1(t, z) + F2(t, z). One part of the vectorfield F1 + F2, say F1, is treated implicitly by the scheme, while the other, i.e., F2, is treated explicitly. Different IMEX schemes have been developed, under which there are the popular IMEX-BDF schemes that are of multistep type [26]. We, however, will use a Runge-Kutta IMEX scheme, because we prefer to work with one-step methods. Especially when working with spatial adaptivity, implementation of multistep methods can become very complicated due to the fact that every time level has its own spatial discretization. We will use an IMEX-midpoint scheme, which can be seen to be a combination of an implicit and an explicit midpoint step, and is given by zs = zn + 12τF1(tn + 1 2τ, zs) + 1 2τF2(tn, zn), zn+1 = 2zs − zn + τ ( F2(tn + 12τ, zs)− F2(tn, zn) ) , (6.13) where the s in zs refers to the intermediate stage. This is a second order time integration method and the implicit part, i.e., the implicit midpoint method, is Astable. Also, using this scheme for the systems at hand never revealed any stability problems. Application of (6.13) In the application of this method to the system (6.2)–(6.6), we use the representation (6.8) and choose the implicit and explicit parts as shown in